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ABSTRACT 

We examine the hypothesis that there exists a simple scaling between the observed 
velocities of jets found in Young Stellar Objects (YSOs) and jets found in Active Galac¬ 
tic Nuclei (AGN). We employ a very simplified physical model of the jet acceleration 
process. We use time-dependent, spherically symmetric wind models in Newtonian 
and relativistic gravitational fields to ask whether the energy input rates required to 
produce the jet velocities observed in YSOs (of about 2 x the escape velocity from 
the central object) can also produce AGN jet velocities (Lorentz factors of 7 ~ 10). 
Such a scaling would be expected if there is a common production mechanism for such 
jets. We demonstrate that such a scaling does exist, provided that the energy input 
process takes place sufficiently deep in the gravitational potential well, enabling phys¬ 
ical use to be made of the speed of light as a limiting velocity, and provided that the 
energy released in the accretion process is imparted to a small fraction of the available 
accreting material. 
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1 INTRODUCTION 


Highly collimated jets are observed in a variety of astrophysical objects. They have been found in quasars, active galactic nuclei 
(AGN), stellar binaries, planetary nebulae, young stellar objects (YSOs) and young pulsars (Ferrari 1998; Livic 1999; Reipurth 
& Bally 2001; Gotthelf 2001). However, despite a large theoretical effort, there is still no agreement as to the mechanisms 


which give rise to the acceleration and to the collimation of these jets ( 

Blandford 

1993, 

2000, 

2002; Pringle 

1993; 


1998; Livio| 1999). In this paper we shall restrict our discussion mainly to jets in AGN and in YSOs, although our findings 


will of course have relevance to other areas. In both AGN and YSOs, inflow of matter is thought to be occurring through 
an accretion disc, and it is this process which is thought to power the jets. Also in these objects, the absence of substantial 
thermal emission implies that the jets are not simple hydrodynamic flows, powered by thermal pressure (discussed for example 
by Blandford & Rees 1974; Konig] 198^ ). For this reason it is usually assumed that the acceleration mechanism is associated 
with magnetic fields. Because the first collimated jets to be discovered were the relativistic radio jets from galactic nuclei, 


presumably powered by accretion onto the central black hole (Rees 1984), all the early models of jets involved the generation 


of jets from a central black hole (Ferrari 1998). Indeed the frequently invoked Blandford-Znajek mechanism (Blandford & 


Znajek 1977) for such jet formation requires the tapping of rotational energy from a rapidly spinning black hole. It is clear, 


however, that such exotic processes are of little relevance to the generation of jets around young stellar objects. This leaves 
jet theorists with a dilemma. Do we argue that the jet generation mechanisms in AGN and YSOs are unrelated, and that all 
the various jet formation models put forward so far are correct when applied to the right object? Or can we apply Occam’s 
razor, and argue that all jets are produced by fundamentally the same mechanism, and that essentially the same theory can 


be applied in both cases, when the appropriate scalings are applied (Konig) 1986; Pringli 


1993 


Livio 


19971, |199S|)? 


The jet acceleration process by a spinning disc is often envisaged as being due to centrifugal acceleration of disc material by 


large scale poloidal field lines threading the disc ( 

Blandford & Payne 

1982, Pudritz & Norman 

1986). This idea is exemplified 

by the magnetic wind solution of 

Blandford & Payne 1 

1982) 

and has been demonstrated in 

a number of numerical simulations 

(Ouyed & Pudritz 

1997 

1999 

; Duyed, Pudritz & Stone 

1997 

Kudoh, Matsumoto & Shibata 

1998; Koide et al. 

200c|). It appears 


to be a fairly generic observed property of collimated jets that the jet velocities are comparable to the escape velocities from 























































































































2 Price, Pringle and King 


the central gravitating objects. In the YSO jets, the intrinsic jet velocity is hard to measure, because in order to see the 
jet it needs to have interacted with some of the surrounding material, and so to have been slowed to some extent. Thus 
it is difficult to measure the core of the jet directly. Typically the jet velocities at various part of the flow are inferred by 
modelling the velocity structures in the neighbourhood of HH objects (which are basically shocks within the jet) and from 


the observed proper motions of the HH objects, which of course give lower limits to the jet velocity ( 
Such considerations indicate that typical YSO jet velocities are in the range Vjet ~ 300 — 500 krn/s ( 

Leipurth & Bally 

2001 ). 

Eisloffel & Mundt 

1998; 

Micono et al. 1998; Bally et al. 2001; Hartigan et al. 2001; Reipurth et al. 2002; Bally et al. 2002; 

Pyo et al. 2002) con 

ipared 

km/s. 

to the escape velocity from a typical young star (mass 1 Mg, radius 5 Rq; Tout, Livio & Bonne! 

19990 ofiw ~ 270 


The AGN jets are relativistic, appropriate for a velocity of escape from close to a black hole, and appear observationally to 
have relativistic gamma-factors around 7 j et ~ 5 — 10 (Urry & Padovani 1995; Biretta, Sparks & Macchettc 199£), although 
arguments for higher values ( 7 j e t ~ 10 — 20) have been made on theoretical grounds (Ghisellini & Celott: 2001). It is clear 
from this that scale-free models for driving outflows from discs, such as the model of Blandford & Payne (1982), are lacking 
in the basic respect that the observed disc-launched jets appear to know of the scale associated with, and thus presumably to 
be launched from, their very central regions ( Konigj 1986; Pringle 1993; Livic 1997). In addition, it now seems unlikely that 


accretion discs are able to drag in the large-scale poloidal fields assumed in the models to thread the inner disc regions (Lubow. 


Papal hzou & Pringle 1994). This has led to a number of recent suggestions that jet acceleration, and perhaps collimation, 


can be due to locally generated, small scale, perhaps tangled, magnetic fields ((Tout & Pringle ! 

996; 

Turner, Bodenheimer & 

Rozyc 

zka 

1999; 

Heinz & Begelman 

200 C; 

Kudoh, Matsumoto & Shibata 

2002 ; L 

2002 ; 

William 

s 2002 ). 


Bearing this in mind, we set out here to address one specific question, which is: are the acceleration mechanisms for 
the jets the same in both AGN and YSOs? Given the uncertainties of the jet acceleration process itself, we approach the 
problem in a rather generic and abstract manner. Since we are only concerned here with the acceleration process, and not the 
jet collimation, we consider the driving of a spherically symmetric outflow by the injection of energy into the gas at a fixed 
radius close to the central object^]. This, in general terms, must be the basis of any jet acceleration mechanism. And, since 
the final jet velocities are directly related to the size of the central object, it follows that the acceleration process must be 
reasonably well localised in that vicinity. We treat the gas in a simple manner as having a purely thermal pressure, P, and 
internal energy, it, and a ratio of specific heats 7 which we shall take to be 7 = 4/3. The exact value of 7 is not critical to 
our arguments, provided that 7 < 5/3 so that the outflow becomes supersonic. We note, however, that taking 7 = 4/3 is in 
fact appropriate to the case of an optically thick radiation-pressure dominated flow, and to the case in which the dominant 
pressure within the gas is caused by a tangled magnetic field (e.g Heinz & Begelman 2000). Thus we feel that such treatment 
should allow us to draw some general conclusions. 

If the jet acceleration mechanism is the same for both YSOs and AGN, then the same (appropriately scaled) energy input 
rate should account for the jet velocities in both sets of objects. Thus, for example, we ask whether the same energy input 
rate gives rise to a final jet velocity Vj e t ~ 2v eS c in the non-relativistic case and 7j e t ~ 7 in the relativistic case. With this 
in mind, we undertake the following computations. In Section 2, we consider non-relativistic outflows, relevant to YSOs. In 
Section 3, we consider relativistic outflows, from compact relativistic objects, relevant to AGN. In both cases we start from a 
steady configuration, in hydrostatic equilibrium, and inject energy into the gas at a steady rate over a small volume close to 
the inner radius. We follow the time-evolution of the gas as it expands. Once the expansion has proceeded to a large enough 
radius, we match the solution onto a steady wind solution in order to estimate the final outflow velocity. We then plot the final 
jet velocity as a function of the (dimensionless) energy input rate (heating rate) for both the relativistic and non-relativistic 
cases. 

In Section U, we present our results and conclusions. 


2 NON-RELATIVISTIC (YSO) JETS 


2.1 Fluid equations 


For YSO jets we expect the gravitational field to be well approximated by a non-relativistic (Newtonian) description. In one 
(radial) dimension the equations describing such a fluid including the effects of energy input are expressed by the conservation 
of mass, 


dp r dp . P d 2 rs n 
+ v — + —^ r v ) = 0 ’ 


dt 


dr 


momentum, 


(1) 


1 Note that we could equivalently consider an injection of momentum rather than energy. This might add physical reality at the expense 
of complexity, and here we choose to remain with our simplified abstract approach 
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dv r r dv r 1 dP GM 

dt + dr p Or r 2 

and energy, 


d(pu) r d(pu) 
dt V dr 


+ 


P + pu 


^(rV) = pA, 


( 2 ) 

(3) 


where p, v r , P and u are the fluid density, radial velocity, pressure and internal energy per unit mass respectively, M is the 
mass of the gravitating object (in this case the central star), and 


dQ _ T <ls_ 
dt dt 


(4) 


is the heat energy input per unit mass per unit time (where T and s are the temperature and specific entropy respectively). 
The equation set is closed by the equation of state for a perfect gas in the form 


P = {7 - !)p u - 


(5) 


2.1.1 Scaling 

To solve ( 0 )-® numerically we scale the variables in terms of a typical length, mass and timescale. These we choose to be 
the inner radius of the gas reservoir [L] = R f , the mass of the gravitating body [M] = M * and the dynamical timescale at 
the the inner radius (r = R*), [t] = (GM,/i^)"^ 2 . In these units GM = 1 and the density, pressure, velocity and internal 
energy, respectively, have units of density, [p] = M,/Rl, pressure, [P] = M*/(P»r 2 ), circular velocity at f?.*, [n] = \J GM,/R * 
and gravitational potential at R t , [ m ] = Note that the net heating rate per unit mass A is therefore given in units 

of gravitational potential energy, GM,/R *, per dynamical timescale at R*, (GM*/P*) -1 ^ 2 . We point out that this scaling 
is simply to ensure that the numerical solution is of order unity and that when comparing the results to the relativistic 
simulations we scale the solution in terms of dimensionless variables. 


2.2 Numerical solution 


We solve ( 0 )-(§) in a physically intuitive way using a staggered grid where the fluid velocity is defined on the half grid points 
whereas the density, pressure, internal energy and heating rate are specified on the integer points. This allows for physically 
appropriate boundary conditions and allows us to treat the different terms in a physical way by applying upwind differencing 
to the advective terms but using centred differencing on the gradient terms. The scheme is summarised in Figure A1 with 
the discretized form of the equations given in appendix |a| The staggered grid means that only three boundary conditions 
are required, as shown in Figure Al. We set v = 0 at the inner boundary and the density and internal energy equal to their 
initial values (effectively zero) at the outer boundary. 


2.3 Initial conditions 


The form of the initial conditions is not particularly crucial to the problem, as the wind eventually reaches a quasi-steady 
state that is independent of the initial setup. What the initial conditions do affect is the time taken to reach this steady state 
(by determining how much mass must initially be heated in the wind). We proceed by setting up a body of gas (loosely ‘an 
atmosphere’) above the ‘star’ (or rather, an unspecified source of gravity) initially in hydrostatic equilibrium, such that v = 0 
everywhere and 


dP_ _ GMp 
dr r 2 

The pressure is related to the density by a polytropic equation of state 


( 6 ) 


P = Kp\ 


(7) 


where K is some constant. Combining these two conditions we obtain an equation for the density gradient as a function of 
radius 


dp(r) _ p{r)-^~ 2) GM 

dr 'yK r 2 


Integrating this equation from r to some upper bound R a 


p{r) = 


7-1 (GM GM\ 

I- ?TV“ 7W 


i/(7-i) 


we obtain 


(9) 
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[h] 



1 10 100 1000 


Figure 1 . The initial conditions for the non-relativistic case, We plot profiles of density, pressure and intern al ener gy per unit mass (or 
temperature) as functions of radius. The quantities here are dimensionless and the units are as described in §^.l.l]. 


To ensure that pressure and density are finite everywhere (for numerical stability) we set Roo = oo. The density is then given 
as a simple function of radius where it remains to specify the poly tropic constant K. In scaled units we choose K = ( 7 — l )/7 
such that p(R „) = 1 ( i.e. the central density equals the mean density of the gravitating body - note that we neglect the 
self-gravity of the gas itself. Choosing 1\ effectively determines the amount of mass present in the atmosphere and thus the 
strength of the shock front which propagates into the ambient medium (in terms of how much mass is swept up by this front). 

We set the initial pressure distribution using (Q). If we do this, however, the slight numerical imbalance of pressure and 
gravity results in a small spurious response in the initial conditions if we evolve the equations with zero heating. In the 
non-relativistic case the spurious velocity is kept to an acceptably small level by the use of a logarithmic radial grid (thus 
increasing the resolution in the inner regions). In the relativistic case however this slight departure from numerical hydrostatic 
equilibrium is more significant. This response is therefore eliminated by solving for the pressure gradient numerically using the 
same differencing that is contained in the evolution scheme. That is we solve from the outer boundary condition P(r ma x) = 
Kp(r m ax ) 7 according to 


Pi -1 = Pi - (r; - n -1 


Pi-1/2 


' i- 1/2 


( 10 ) 


Solving for the pressure in this manner reduces any spurious response in the initial conditions to below round-off error. The 
internal energy is then given from (Q). The pressure calculated using ([to]) is essentially indistinguishable from that found 
using (,^)(A P/P ~ 1CU 5 ). The initial conditions calculated using equation ( 1 ), © and (joj) are shown in Figure [I]. We use 
a logarithmic grid with 1001 radial grid points, setting the outer boundary at r/R » = 10 3 . Using a higher spatial resolution 
does not affect the simulation results. 


2.3.1 Heating profile 

The choice of the shape of the heating profile A(r) is fairly arbitrary since we wish simply to make a comparison between 
the non-relativistic and relativistic results. We choose to heat the wind in a spherical shell of a fixed width using a linearly 
increasing and then decreasing heating rate, symmetric about some heating radius rh ea t which we place at r = 2.IT?*. The 
heating profile is spread over a radial zone of width 2 R* (that is the heating zone extends from r = 1.17?* to 3.1i?*)(see 
Figure ^|). We choose a heating profile of this form such that it is narrow enough to be associated with a particular radius of 
heating (necessary since we are looking for scaling laws) whilst being wide enough to avoid the need for high spatial resolution 
or complicated algorithms (necessary if the heat input zone is too narrow). The important parameter is thus the location of 
the heating with respect to the Schwarzschild radius, so long as the heating profiles are the same in both the relativistic and 
non-relativistic cases. Provided that the heating profile is narrow enough to be associated with a particular radius and wide 
enough to avoid numerical problems, our results do not depend on the actual shape of the profile we choose. 
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1 10 100 1000 


radius 

Figure 2. Results of a typical non-relativistic simulation at time t = 1000 (where units of time are the dynamical time at the innermost 
radius, yjR^/GM). Quantities shown are the Mach number (y/c a ), velocity, heating rate (A), internal energy per unit mass (u = rttherm), 
log(density) and log(pressure). 


2.4 Results 

The results of a typical non-relativistic simulation with a moderate heating rate are shown in Figure [2] at £ = 1001 (where 
t has units of the dynamical time at the inner radius). We observe the effect of the heating propagating outwards in the 
atmosphere in the form of a shock front. After several hundred dynamical times the wind structure approaches a steady state 
in that there is only a small change of the overall wind structure due to the shock continuing to propagate outwards into 
the surrounding medium. The small disturbance propagating well ahead of the main shock is a transient resulting from the 
response of the atmosphere to the instantaneous switch-on of the heating. The velocity of the gas begins to asymptote to 
a constant value as the shock propagates outwards. Plotting the mass outflow rate M = 4nr 2 pv and the Bernoulli energy 
E = ip 2 + pu + P — GM/r as a function of radius (Figure |d|), we see that indeed the wind structure is eventually close to 
that of a steady wind above the heating zone (ie. M and E ~ constant). It is thus computationally inefficient and impractical 
to compute the time-dependent solution for long enough to determine an accurate velocity as r —> 00 when the wind will 
continue to have a steady structure. Instead we find the steady wind solution for a given amount of energy input to the wind 
corresponding to the energy plotted in Figure^ (top panel). 


2.5 Steady wind solution 


Non-relativistic, steady state ( d/ dt = 0) winds with energy in put have been well studied by ma ny authors, and t h e equ ations 
describing them can be found in Laniers fc Cassinelli (1999), who credit the original work to Holzer fe Axford (197()|). The 


reader is thus referred to 


Lamers & Cassinelli 


( 1999 ) for details of the derivation. As in the usual Bondi/Parker ( Bondi 1955, 
Parker 1958 ) wind solution with no heat input, we set d/dt = 0 in ( 0 )-® and combine these equations into one equation for 
the Mach number M 2 = v 2 /c 2 as a function of radius, given by 
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1 10 100 1000 

radius 

Figure 3. Bernoulli energy E = ^v 2 + pu + P — GM/r (top) and mass outflow rate M = 4nr 2 pv (bottom) in the time-dependent wind 
solution at t = 1000. The profiles are approximately constant over the region between the two circles. The sample point used to match 
this flow to the appropriate steady state solution is indicated by a cross. 


W 



1 10 100 1000 10 4 


radius 

Figure 4. Steady wind Mach number (top panel) and velocity (centre panel) profiles are compared to the time-dependent solution (plotted 
every 100 dynamical times). There is a small discrepancy between the two solutions where we have taken the limit in approaching the 
singular point at M = 1, but an otherwise excellent agreement between the two solutions. 


dM 2 = Af 2 (2 + ( 7 ~1)A/ 2 ) [ 2 dQ GM (5 - 3 7 ) _ 4e(r) 

dr 2(Af 2 — 1) [e(r) + GM/r] ^ dr r 2 ( 7 — 1) r 


( 11 ) 


where dQ/dr is the local heating gradient and e(r) is the Bernoulli energy which is specified by integrating the Bernoulli 
equation 


de(r) 

dr 


d 

dr 


1 2, , p GM 

-v +pu + P - 

2 r 


dQ 
dr ’ 


to give 
e(r) = 


e(roo) - Q(r) 


r 


( 12 ) 




















dQ 
dr ’ 
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(13) 


where Q(r) is the total energy input to the wind. Since we are interested in the terminal velocity of the outflow we choose 
a point above the heating shell where the energy has reached its steady state value (i.e. where the energy is constant in 
Figure [lj top panel) and integrate outwards using the energy and Mach number at this point to solve (|ll|) as an initial value 
problem. Note that in fact the terminal velocity is determined by the (constant) value of the Bernoulli energy above the 
heating zone since as r —> oo, e(r) —> ^v 2 . However we compute the steady wind profiles both inwards and outwards to show 
the consistency between the time-dependent solution and the steady state version. 

In order to perform the inward integration, we must determine the energy at every point for our steady solution by 
subtracting the heat input from the steady state energy as we integrate inwards through the heating shell (^) . To determine 
this however we must also determine the local (steady state) heating gradient dQ/dr, which is related to the (time dependent) 
heating rate A by setting d/dt = 0 in the time dependent version, ie. 


. dQ dQ dQ dQ 
A = -f- = +v-A = v-?-. 
dt dt dr dr 

We therefore calculate dQ/dr from the time dependent solution using 


(14) 


dQ _ A(r) 
dr v(r) ’ 


(15) 


where v(r) is the wind velocity at each point in the heating shell from the time-dependent solution. The problem with this is 
that at the inner edge of the heating shell the heating rate is finite while the velocity is very close to zero, resulting in a slight 
overestimate of the total energy input near the inner edge of the shell in the steady wind solution. Care must also be taken 
in integrating through the singular point in equation ( |lT[ ) at M 2 = 1. Most authors (e.g. Hamers & Cassinelli 1999) solve the 
steady wind equations starting from this point but for our purposes it is better to start the integration outside of the heating 
shell where the energy is well determined. We integrate through the critical point by using a first order Taylor expansion and 
appropriate limit(s), although this introduces a small discrepancy between the steady state and time-dependent results in this 
region (Figure^). 

Having determined the energy and heating gradient at each point in the wind we integrate ( [TT] ) both inwards and outwards 
from the chosen point above the heating shell using a fourth order Runge-Kutta integrator (scaling (Jn|) to the units described 
in §|2.1.l[). The velocity profile is then given by v 2 = M 2 c 2 where 


l(r) = 


2(7 ~ 1) 


2 + M 2 (r)(y — 1) 


e(r) + 


GM 


(16) 


The resulting steady wind solution is shown in Figure [f] along with the time-dependent solution. The two profiles are in 
excellent agreement, proving the validity of our time-dependent numerical solution and the assumption that the wind is in 
a steady state. The steady solution thus provides an accurate estimate of the velocity at arbitrarily large radii (although as 
pointed out previously this is set by the value of the steady state Bernoulli energy). 


2.6 Terminal wind velocities as a function of heating rate 

Using the steady wind extrapolation of the time-dependent solution, we can determine the relationship between the heating 
rate and the terminal wind velocities. In order to make a useful comparison between the heating rates used in both the 
Newtonian and the relativistic regimes, we need to define a local canonical heating rate A c (r) valid in both sets of regimes. In 
dimensional terms the heating rate A(r) corresponds to an input energy per unit mass per unit time. Thus we need to define 
the local canonical heating rate as 

A E 

A c (r) = —, (17) 

for some relevant energy A E and some relevant timescale At. 

There are clearly many different ways in which we might define a canonical heating rate. We find, however, that our results 
are not sensitive to the particular choice we make. We shall make use of a definition which draws on the physical processes 
we expect to be behind the jet acceleration process. Even though the processes by which this occurs are still obscure, we 
expect the energy for the jet to be provided fundamentally by liberation of energy in a rotating flow. Thus, with this physical 
motivation in mind, we take the canonical energy per unit mass, A E, to be the energy released locally by bringing to rest 
a particle of unit mass which is orbiting in a circular orbit at radius r. In the Newtonian regime this is simply the kinetic 
energy of a circular orbit 
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[t] 



Average dimensionless heating rate < A > 

Figure 5. Terminal wind velocities plotted as a function of the average dimensionless heating rate (A). Wind velocities are plotted in 
units of the escape velocity at the inner radius (ie. r = ft* = 1), Vesc = (2GM/ft*) 1 / 2 . We compute solutions corresponding to velocities 
typically observed in YSO jets (with a fairly generous upper limit of v/v esc ~ 3). 


(An alternative possibility, for example, would be to take A E to be the energy released by dropping a particle from infinity and 
bringing it to rest at radius r, which would correspond to the escape energy from that radius, GM/r.) By similar reasoning, 
we take the canonical timescale on which the energy is released to be the orbital timescale at radius r, that is At = fi^ 1 , 
where 

Q 0 = (GM/r 3 ) 1/2 . (19) 


Using this, we are now in a position to define a local canonical heating rate as 

A c (r) = A £xll„=^. (20) 

For intercomparison of our various wind computations both in the Newtonian and in the relativistic regimes, we now use the 
canonical heating rate derived above to define a dimensionless heating rate for each wind computation. Because heat is added 
over a range of radii, we need to define the dimensionless heating rate (A) as an appropriate volume average. We shall define 


(A) 


f r2 A r 2 dr 

J r\ 

A c (r max ) r 2 dr ’ 


( 21 ) 


where r ma x is the radius at which the heating rate A(r) takes its maximum value and ri and V 2 are the lower and upper 
bounds of the heating shell respectively. 

The relation between this average dimensionless heating rate and the terminal wind velocity is shown in Figure The 
wind velocities are plotted in units of the escape velocity v e sc at R, and solutions are computed for wind velocities of up to 
~ 3v esc - The important point in the present analysis is that the heating rate can be meaningfully compared to the relativistic 
results (see below). 


3 RELATIVISTIC JETS 

Having determined the heating rates required to produce the observed velocities in YSO jets we wish to perform exactly the 
same calculation within a relativistic framework. We proceed in precisely the same manner as in the non-relativistic case. We 
adopt the usual convention that Greek indices run over the four dimensions 0,1,2,3 while Latin indices run over the three 
spatial dimensions 1,2,3. Repeated indices imply a summation and a semicolon refers to the covariant derivative. The density 
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p refers to the rest mass density only, that is p = nmo where n is the number density of baryons and mo is the mass per 
baryon. 


3.1 Fluid equations 

The equations describing a relativistic fluid are derived from the conservation of baryon number, 

(PUA;» = 0 , 


( 22 ) 


the conservation of energy-momentum projected along a direction perpendicular to the four velocity (which gives the 
equation of motion), 


h„ a T a % = (g„ a + U^U a )T a G = 0 , 

and projected in the direction of the four-velocity (which gives the energy equation), 

U a T a G = 0 . 

Here the quantity T A “' is the energy momentum tensor, which for a perfect fluid is given by 

c 2 T^ v = phU^U 1 ' + Pg^, 

where h is the specific enthalpy, 

J. 2 , p 2 , 7 P 

h — c + m + — — c + --r—. 

p (7 - 1 )p 

As in the non-relativistic case u is the internal energy per unit mass, P is the gas pressure and we have used the equation of 
state given by equation (jjij). The energy equation may also be derived from the first law of thermodynamics using equation (p2|), 
which is a more convenient way of deriving an energy equation in terms of the internal energy (rather than the total energy) 
and in this case ensures that the meaning of the heating term is clear. The metric tensor is given by the Schwarzschild 
(exterior) solution to Einstein’s equations, that is 

-l 


(23) 

(24) 

(25) 

(26) 


ds 2 = -A dr = - 1 - 


2 GM \ 
c 2 r ) 


c dt + 1 - 


2 GM Y 
c 2 r ) 


dr + r(d9 2 + sin 6d(j> )• 


We consider radial flow such that U° = = 0. The four velocity is normalised such that 


U„l = -c 2 , 
and we define 

U* = = ( 1 - 

dr 


2GM 


c 2 r J 


1 - 


2 GM\ ( U r ) 2 


I 


1/2 


which we denote as 


^ = 4 


where we set for convenience 


r = 

and 


1 - 


2 GM\ x (U r ) 
c 2 r ) 


21 1/2 


+ ' 


a = 1 - 


2GM 


■ ) 
c z r J 


(27) 


(28) 


(29) 


(30) 


(31) 


(32) 


Note that while a corresponds to the lapse function in the 3+1 formulation of general relativity, the quantity T is not the 
Lorentz factor of the gas (which we denote as W) as it is usually defined in numerical relativity (e.g. Banyuls et al. 1997) but 
is related to it by W = T /a. From (|2f]) we also have the relation 

dUt ^ ^ (33) 


dt a 2 Tc 2 dt 

From (^), (^3|) and (|2^) using (^), (^]), (^) and (|3^) we thus derive the continuity equation, 

U r dU r ~ 


dp T dp a 2 p 

i +v i + -T- 


— — < r 2 U r ) + 
r 2 dr' a 2 Tc 2 dt 


= 0, 


(34) 


the equation of motion, 
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dU 


r dU r ro?c 2 dP ir_dP_ a 2 GM _ 

dt ~ >rV dr ^ ph dr ph dt F r 2 ’ 


and the internal energy equation, 
d(pu) r d(pu) a 


dt 

where 


+ v 


dr 


+ — (P + pu) 


1_ d_ 

r 2 dr 


(r 2 U r ) + 


U r dU r 


a 2 rc 2 dt 


a A 

= — pA, 


U r 




_ dr 

U t= Yt 

is the velocity in the co-ordinate basis. We define the heating rate per unit mass as 
, ds 


A = T 


dr ’ 


(35) 


(36) 


(37) 


(38) 


where T is the temperature, s is the specific entropy and dr refers to the local proper time interval (A is therefore a local rate 
of energy input, caused by local physics). A comparison of (|3j|), ( ^B| ) and (^) with their non-relativistic counterparts (jj|), (|2j) 
and (j|) shows that they reduce to the non-relativistic expressions in the limit as c —* oo, and to special relativity as M —> 0. 

The ‘source terms’ containing time derivatives of U r and P are then eliminated between the three equations using the 
equation of state (|E]) to relate pressure and internal energy. Substituting for pressure in ( p6| ) and substituting this into ( |35| ) 
we obtain the equation of motion in terms of known variables, 


dU 


dt X \ ph 


dU r _ _ 

dr phTX dr 

where for convenience we define 
'7 P\ U r U r 

ph J r 2 c 2 ’ 


GM v r 7 P 2 U r 


rx 


+ X ph r hX ^ A ’ 


X = 1 — 


(39) 


(40) 


and we have expanded the \j§z(r 2 U r ) terms in order to combine the spatial derivatives of U r into one term. We then 
substitute (Eh}) into (E3) and (liq) to obtain equations for the density 


dp r dp or 
dt +V d^ ~ 


pA- 


dP U r U r (7-1) 


hVX dr r 2 c 2 X 


pA 


and internal energy, 


d(pu) 


+ v r 1- 


7 P a 2 \ d(pu) 


dt V ph T 2 A ') dr 

where for convenience we have defined 


a 

Y 


(P + pu)A — 1 + 


U r U r 7 P 
T 2 c 2 X ph 


pA 


A = 


1 - 


U r U r 

r 2 c 2 x 


1 - 


rP 

ph 


dU r 

dr 


+ 


U r U r 

+ rw 1 ~ 


rP 

ph 


2 U r 
r 


U r GM 


(41) 


(42) 


(43) 


r 2 c 2 A r 2 

From the solution specifying U r we calculate the velocity measured by an observer at rest with respect to the time slice 
(referred to as Eulerian observers), which is given by 

U r v r 

* r = YP= a ’ (44) 

since there are no off-diagonal terms (ie. zero shift vector) in the Schwarzschild solution. For these observers the Lorentz factor 
is given by 

-1/2 


W= 1 - 


(45) 


where v r v r = g rr v r v r , such that U r = Wv r . 


3.2 Scaling 


The usual practice in numerical relativity is to scale in so-called geometric units such that G = M = c = 1. In these units the 
length scale would be the geometric radius GM/c 2 and the velocity would have units of c. Instead for the current problem, 
we adopt a scaling analogous to that of the non-relativistic case, that is we choose the length scale to be the radius of the 
central object, R*, where R „ is given as some multiple of the geometric radius, ie. 


m p GM * 

[L\ =R*= n—— 


with n ^ 2.0. The mass scale is again the central object mass [M] = M* and the timescale is given by 


(46) 
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M 



Figure 6. The initial conditions for the gas reservoir for the relativistic cases of a neutron star (dashed line) ( R*/Rs c h = 2.5) and 
black hole (solid line)(i?*/i?Sch = 1.0). Note, however, that the innermost radius is at r = l.Oli?* in the latter case. We plot profiles 
of internal energy per unit mass (or temperature), density and pressure, as functions of radius. These quantities are given in units of 
GM/i?*, M/R * and M* / (R*t1) respectively. Note that steeper gradients are required to hold the gas in hydrostatic equilibrium as the 
gravitational field becomes more relativistic. The black hole reservoir is of lower density than the neutron star version because of the 
choice of the polytropic constant (chosen such that the central density is of order unity). 


GM.V 172 3/ 2 GM, 

m ) ~ n c 3 


(47) 


In these units, velocity is measured in units of [u] = n _1,/2 c (or equivalently c 2 = n). The scaled equations are thus given 
simply by setting G = M = 1 and c 2 = n everywhere. 

This scaling ensures that the relativistic terms tend to zero when c (or n) is large and that the numerical values of p, 
pu and U r are of order unity. We thus specify the degree to which the gravity/gas dynamics is relativistic by specifying the 
value of n (i.e. the proximity of the innermost radius, and thus the heating, to the Schwarzschild radius, i?sch = 2 GM/c 2 ). 
We compute solutions corresponding to gas very close to a black hole (highly relativistic, n = 2.0, or R, = 7?s c h), neutron 
star (moderately relativistic, n = 5, or R, = = 5 GM/c 2 , which is equivalent to heating further out and over a wider 

region around a black hole) and white dwarf/non-relativistic star (essentially non-relativistic, n = 5000, or R, = 2500-Rsch)- 
Note that in the highly relativistic case although we scale the solution to n = 2.0 such that the mass, length and time scales 
(and therefore the units of heating rate, energy etc.) correspond to those at r = 7?sch, our numerical grid cannot begin at R * 
as it does in the other cases. We therefore set the lower bound on the radial grid to slightly below the heating shell (typically 
r = 1.01.R* where the heating begins at 1.11?*). Note that the above scaling is merely to ensure that the numerical solution 
is of order unity, since we scale in terms of dimensionless variables to compare with the non-relativistic solution. 


3.3 Numerical Solution 


In order to s olve the relativistic fluid equations numerically we use a method analogous to that used in the non-relativistic 
case (Figure Al). That is, we first compute \J T on the staggered (half) grid and use this to solve for p and pu on the integer 
grid points. Again the advective terms are discretized using upwind differences (where the ‘upwindedness’ is determined from 
the sign of the co-ordinate velocity v r ) and other derivatives are calculated using centred differences. As in the non-relativistic 
case, where a centred difference is used, the quantities multiplying the derivative are interpolated onto the half grid points if 
necessary. In equation (hit) we evaluate the 8P/dr term using upwind differences. 


3.4 Initial Conditions 

We determine initial conditions for the relativistic case by setting U r = 0 and d/dt = 0 in ([fi]), from which we have 

2 GM 
c 2 r 


dP_ 

dr 


ph GM 


1 - 



(48) 
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The pressure is thus calculated as a function of p, u and P (where P = (7 — 1 )pu). We solve ( [f 8 | ) using the same assumptions 
as in the non-relativistic case (§[ 2 ^), that is an adiabatic atmosphere such that 

P = Kp 7 . (49) 


We therefore have 

'K = _L_ + 7 K P 

dr ■yKa 2 c 2 ( 7 —1) 


GM 
r 2 ’ 


(50) 


which we solve using a first order (Euler) discretization to obtain a density profile. The pressure may then be calculated using 
([tol), however to ensure that hydrostatic equilibrium is enforced numerically we solve ([fi| ) using the same discretization as in 
the fluid equations, integrating inwards from the outer boundary condition P(r max ) = A'p(r ma x) 7 - However in this case the 
pressure gradient also depends on the pressure, so we use the pressure calculated from ( (til) to calculate the initial value of the 
specific enthalpy h and iterate the solution until converged ([ P n+1 — P n yp n < 10~ 10 ). In the black hole case the resulting 
pressure differs from that found using (^) by A P/P ~ ICC 2 . We choose K such that the central density is of order unity - 
typically we use K = 107/(7 — 1) in the black hole case. Note that changing K simply changes the amount of matter present 
in the atmosphere but does not affect the temperature scaling and does not affect the final results (although it significantly 
affects the integration time since it determines the strength of the shock front and the amount of mass to be accelerated). 

Initial conditions calculated in this manner for the black hole ( R^/Rsch = n/2 = 1.0) and neutron star (R„/Rsch = 2.5) 
atmospheres are shown in Figure []. The initial setup reduces to that of Figure |l| in the non-relativistic limit when the same 
value of K is used. We set the outer boundary at r/R » = 10 4 , using 1335 radial grid points (again on a logarithmic grid). 


3.5 Results 


The results of a typical (n=2.0) relativistic simulation are shown in Figure ^ at t = 1000. Again we observe that the wind 
structure reaches a quasi-steady state, with the velocity approaching a steady value at large radii. Note that because the 
steady state density is higher than that of the surrounding medium no wide shock front is observed. 

Plotting the mass outflow rate M = Airr 2 pU r and the relativistic Bernoulli energy E le 1 = ^ V 2 h 2 /c 2 — ^c 2 (see e.g. 


Shapiro & Teukolsky 1982) as a function of radius (Figure | 8 |), we see that indeed the structure approaches that of a steady 
(relativistic) wind (that is, the energy and M profiles are flat above the heating zone). We may thus apply a relativistic steady 
wind solution with this Bernoulli energy as an initial value to determine the final velocity and Lorentz factor as r —> 00 . 
Note that we cannot apply a non-relativistic steady wind solution because although the gravity is non-relativistic, the outflow 
velocities are not. As in the non-relativistic case the final wind velocity is determined by the steady Bernoulli energy, since in 
this case as r —> 00 , E le 1 —> \[(U r ) 2 — c 2 ]. 


3.6 Steady wind solution 


Relativistic, steady state ( d/dt = 0) winds were first studied by Michel (1972) and extended to include energy deposition by 
Flammang (1982). The problem has recently received attention in the context of neutrino-driven winds in gamma-ray burst 
models by Pruet, Fuller fe Cardalj (2001) and Thompson, Burrows & Meyer (2001). We proceed in a manner analogous to 
that of the non-relativistic solution. Setting d/dt = 0 the continuity (^2[) and momentum ( |23| ) equations become 

1 dU r 


}_dp 

p dr U r 


dr 


2 

4— 
r 


= 0 


= 0 


du r r 2 c 2 dp gm 

dr + ph dr + r 2 
where (|F^) is equivalent to 
r 2 pU r = const. 

Combining (^2|) and (§) we obtain 


1 

TP 


2ti2 2 

r 2 _ C s 

h'y 


dU r 

dr 


c 2 P 2 dc 2 


c 2 r 2 2c 2 


GM 


h'y dr h'y r r 2 ’ 

where c 2 = 7 P/p and (U r ) 2 = U r U r . From the first law of thermodynamics and (| 
equation in the form 


d 

dr 


1 V 2 h 2 

2~P~ 


hr 2 dQ 
c 2 dr ’ 


(51) 

(52) 

(53) 

(54) 

we derive the relativistic Bernoulli 

(55) 


such that both sides reduce to their non-relativistic expressions as c 
as in the non-relativistic case. Expanding this equation we find 


00 . The quantity dQ/dr is the local heating gradient 
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Figure 7. Results of a typical black hole relativistic simulation at t=1000 (where units of time are the dynamical time at the central 
object). Quantities shown are the Mach number (u/c s ), velocity for Eulerian observers (u r ), heating rate (A), internal energy per unit 
mass (u = Utherm)j log(density) and log(pressure). Units of velocity are such that c = J2 and as in the non-relativistic case energy has 
units of GM/R *. 


dc 2 


dQ 


h d f l^ r ^ 2 \ _ h GM 


dr c 2 r 2 dr \2 y J c 2 T 2 r 2 
Combining (^i|) and (^j|) and manipulating terms, we obtain an equation for (U r ) 2 , 


2(U r ) 2 


r c 2 r 2 2d , ,.c 2 r ( r dQ\ GM 

-(7-1)— ( r — ) - — 


— ( U r ) 2 = - 

dr K ’ \{U r ) 2 -c 2 r 2 c 2 /h] [ h r 


dr 


(56) 


(57) 


where c 2 and h = c 2 + c 2 /(q — 1) are given functions of known variables by integration of the Bernoulli equation (|5^), in the 
form 


d (v n - r d ® 
dr ^ F dr’ 

to ensure that h does not appear in the heating term on the right hand side. The integration is then 


e(r) = Th = e(roo) — J 


and hence 

h=^l, 


c 2 = (7 - l)(/i - e 2 ). 


(58) 


(59) 


(60) 


The ‘heating gradient’, TdQ/dr, is calculated from the time-dependent solution using 
^dQ t , a(r)A(r) 






















14 Price, Pringle and King 
W 



10 100 1000 10 4 


radius 

Figure 8. The relativistic Bernoulli energy E re \ = ^Th/c 2 — ^c 2 (top) and mass outflow rate M = Anr 2 pU r (bottom) in the time- 
dependent relativistic wind solution with a reasonably high heating rate are shown as functions of radius at time t = 1000. In order to 
match this solution to a steady outflow solution, the Bernoulli energy is assumed to be constant over the region indicated by the two 
circles, and the steady wind solution is computed using initial values at the point indicated by a cross. 


w 



radius 


Figure 9. The radial profiles of the steady wind r-component of four velocity U r (top panel) and of the velocity for Eulerian observers 
v r (centre panel) are compared to the time-dependent solution (plotted every 100 dynamical times) for a typical relativistic calculation 
for the black hole (n=2.0) case. Units are such that c = V2 on the velocity plots. Note the excellent agreement between the two solutions. 


since 


A = T — = < ^- = U* 
dr dr 


dQ r dQ 

~m +v IP 


(62) 


where r is the proper time and U l = Y/a 2 . The velocity profile for an Eulerian observer is then calculated using ( 0 ) and the 
final Lorentz factor Woo using equation (^|). As in the non-relativistic case we choose a starting point for the integration above 
the heating shell and integrate outwards from this point using a fourth order Runge-Kutta integrator in order to determine 
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Figure 10. The terminal r-component of four velocity U r (top panel) and Lorentz factor (bottom panel) of the wind in the non- 
relativistic (o, solid), white dwarf (x, dot- dash ed). neutron star (+, dotted) and black hole (=l dashed) cases, is plotted as a function of 
the dimensionless heating rate defined in §2.6 The top panel may be compared with Figure Win the non-relativistic case. 


the terminal Lorentz factor. The inward integration (and thus the determination of the steady state heating gradient YdQ/dr ) 
is computed only for consistency. We integrate through the singular point in equation (^t]) by taking a low order integration 
with larger steps as this point is approached. 

The solution calculated using is shown in Figure ^[plotted against the evolving time-dependent solution. The profiles 
are in excellent agreement, verifying the accuracy of the relativistic calculation and showing that the wind may indeed be 
described by the steady state solution. 


3.7 Terminal wind velocities and Lorentz factors as a function of heating rate 

In order to compare the relativistic results to those in the Newtonian regime, we define the local canonical heating rate in a 
similar manner to the non-relativistic case, that is 


A c{r) = 


A E 

A t’ 


(63) 
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for some relevant energy A E and some relevant timescale At. As in Section we take the canonical energy per unit mass, 
A E, to be the energy released locally by bringing to rest a particle of unit mass which is orbiting in a circular orbit at radius 
r. For a particle orbiting in the Schwarzschild metric this is the difference, A E, between the energy constants (defined by the 
timelike Killing vector) of a circular geodesic at radius r, and a radial geodesic with zero velocity at radius r. This implies 


(see, for example, Schutz 1990, Chapter 11) 


(64) 


A 77 / 2 1 — 2 GM/rc r or'A/r/ 2 i!/ 2 

AE/C = [1 — “AGM/rc 2 } 1 / 2 - [1 ~ 2GM/r ° 1 ' 

In the Newtonian limit, this reduces to the expected value AE = = GM/'2r. We again take the canonical timescale on 

which the energy is released to be the orbital timescale at radius r as measured by a local stationary observer. For a circular 
geodesic in the Schwarzschild metric, the azimuthal velocity is given in terms of coordinate time, t, by 


d(j>/dt = Q = (GM/r 3 ) 1/2 . 


(65) 


This is the same expression as for the angular velocity of an orbiting particle in the Newtonian limit. But in terms of the 
proper time, r, of a local stationary observer we have, from the metric, 


dr/dt = (1 - 2 GM/rc ) 1/2 , 
and thus d<j>/dr = fl 0 , where 

rl -i 


ni = 


GM 


1 - 


2 GM 


( 66 ) 


(67) 


Using this, the local canonical heating rate is therefore given by 


A c = AE x fi 0 . ( 68 ) 

In the Newtonian limit, r 2 GM/c 2 , this becomes as expected A c — (GM) 3 ^ 2 / 2r 5 ^ 2 . As in the non-relativistic case we use 
the canonical heating rate derived above to define a dimensionless heating rate (A) as an appropriate volume average using 
equation (| 2 l|). 

The final Lorentz factor of the wind plotted as a function of this dimensionless heating rate is given in the bottom panel 
of Figure |l^ in the highly relativistic (black hole), moderately relativistic (neutron star, equivalent to a broader heating shell 
further away from a black hole) and non-relativistic (white dwarf) cases. 

We would also like to make a meaningful comparison of the final wind velocities in units of the escape velocity from 
the star. Note that we cannot simply compare the scaled velocities since we are in effect introducing a ‘speed limit’ in the 
relativistic solution such that the (scaled) relativistic velocity will always be slower than in the equivalent non-relativistic 
solution. However we can compare the velocity for observers along the worldline of a particle in the wind (ie. observers with 
proper time interval dr), U r = dr/dr (that is, the r component of the four velocity, which in special relativity is given by 
U r = 711 7 ', where 7 is the Lorentz factor). Scaling this in units of the (Newtonian) escape velocity from the central object 
(2 GM/R*) 1 / 2 we can make a useful comparison with the non-relativistic results. This velocity is plotted in the top panel of 
Figure [h] against the dimensionless heating rate. 


4 DISCUSSION AND CONCLUSIONS 

We have considered the input of energy at the base of an initially hydrostatic atmosphere as a simple model for the acceleration 
of an outflow or jet. The problem is inherently a time-dependent one, because the flow velocity at the base of the atmosphere 
is zero. Sufficiently large energy input rates give rise to supersonic outflows. We are, of course, unable to compute the outflow 
for an infinite time, and thus cannot directly measure the terminal outflow speed. However, we make use of the fact that if the 
mass in the atmosphere is sufficiently large compared to the mass outflow rate, then at large radii the outflow approximates 
to a steady state (with constant mass flux). We thus match our time-dependent solutions onto steady state outflow solutions 
at large radii and thus determine the terminal velocities for the outflows. We then compute how the terminal velocity of the 
outflow varies as a function of (dimensionless) heating rate in both Newtonian and relativistic gravitational potentials. The 
results are shown in Figures 5 and juj. 

In the top half of Figure [To| we note that dimensionless energy (or momentum) imparted to the outflow at a given value 
of the dimensionless heating rate is larger in the relativistic regime. This comes about simply because a particular element of 
gas cannot escape from the zone in which the heating is occurring with a velocity which exceeds the speed of light. Thus a 
gas element which is accelerated to relativistic energies in the heating zone spends longer in the heating zone than one which 
is not. Indeed from Figure 10, we see that, once the outflows become relativistic, the energy per unit mass in the outflow 
(proportional to 7 ) is proportional to the dimensionless heating rate (A). This comes about because each fluid element spends 
the same time in the heating zone, because it travels through the zone at a velocity ~ c. 
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From Figure 5 we see that a dimensionless heating rate of (A) ~ 17 gives rise to a terminal outflow velocity of Uj e t — 2u e sc 
in a Newtonian potential. From Figure 10, we see that for the same heating rate, the ‘neutron star’ wind, for which the heating 
rate peaks at about 5.2 7?sch becomes mildly relativistic (7j e t ~ 2), whereas the ‘black hole’ wind, for which the heating rate 
peaks at about 2.1 i?sch, leads to an outflow with 7j e t — 11. Similarly a dimensionless heating rate of (A) ~ 55 gives rise to a 
terminal velocity of Uj e t — 3u e sc in the Newtonian case, to an outflow with 7 j e t ~ 4 in the mildly relativistic case, and to an 
outflow with 7jet — 31 in the strongly relativistic case. We have noted above that although the exact numerical values here do 
depend slightly on the exact definition of the dimensionless heating rate, the basic results remain unchanged. For example, 
using the Newtonian dimensionless heating rate (§2.6) in the strongly relativistic case gives a Lorentz factor of 7j e t — 5 for 
the rate which corresponds to Vj e t — 2v eS c in the non-relativistic case. 

We caution that the above analysis does not demonstrate that the model we use provides an adequate description of the 
physics involved in the acceleration process (for example, the jet energy might be initially mainly in electromagnetic form - 
Poynting flux - and only later converted to kinetic energy of baryons) or that there is no intrinsic difference between the jet 
outflows caused by the relativistic nature of the AGN jets (for example, the AGN jets might be mainly in the form of a pair 
plasma and therefore lighter). And we note that it is evident that more detailed physical models need to be developed before 
further conclusions can be drawn. Nevertheless, we suggest that the generic nature of our analysis might give some insight 
into the physical processes involved in the acceleration of jets. 

Thus we conclude, on the basis of the rather simplified physical model we have employed in our analysis, that it is not 
unreasonable to argue that the jets in AGN are simply scaled up, relativistic versions of the jets in YSOs, and that the 
intrinsic jet acceleration mechanism is indeed the same in both the AGN and the YSO contexts. In making this analogy, again 
on the basis of our simplified model for the acceleration process, we find two further physical conditions must hold. First, 
we find that the energy input process, which leads to the acceleration of the outflow, takes place deep in the (relativistic) 
gravitational well for the AGN case (~ 2 Rsch)- What this means physically is that in this way it is possible to make use of 
the limiting velocity, c, to ensure that relativistic fluid elements remain relatively longer in the energy input zone compared 
to their non-relativistic counterparts. Second, we find that the required dimensionless heating rate is much larger than unity. 
The physical implication of this is that the available energy released in the accretion process must be imparted to a small 
fraction of the available accreting material. 
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APPENDIX A: DISCRETIZATION SCHEME FOR NON-RELATIVISTIC EQUATIONS 

The discretization scheme for the non-relativistic fluid equations is summarised in Figure [A l| . Fluxes are calculated on the 
half grid points while the other terms are calculated on the integer points. We solve (|l)-(|5|) in the following manner: The 
numerical equations are solved first for velocity on the half grid points (dropping the superscript r for convenience), 
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where the superscript n refers to the nth timestep and the subscript i refers to ith grid point {v i+1 / 2 ,Pi+ 1/2 thus being 
points on the staggered velocity grid). The quantity pi+ 1/2 is calculated using linear interpolation between the grid points, ie. 
pi+ 1/2 = \{pi + pi+i)- We then solve for the density and internal energy on the integer grid points using the updated velocity, 
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Figure Al. Schematic diagram of numerical method: density and internal energy are defined on the integer points while velocity is 
calculated on the half points. The solution requires one inner boundary condition on v and two outer boundary conditions for p and pu. 
Updated velocities (v n + 1 ) are used to calculate p 71- * -1 and pit 71- * -1 . The scheme allows centred differencing on terms involving staggered 
quantities (top panel) while upwind differencing is used on the advective terms (bottom panel). 
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where At = t n+1 — t n and the timestep is regulated according to the Courant condition 

, min(Ar) 

At < -7T—r+- —T^ 

max(|u|) + max(c s ) 

where c s is the adiabatic sound speed in the gas given by cl = 7 P/p. We typically set At to half of this value. 
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